R Exercises - Particle analysis from SEM images

1 Context

In this exercise we are going to study the results from a particle analysis performed on a set of SEM images of Ni nanoparticles.

To obtain these nanoparticles, we start from Si wafers on which a 5 or 10 nm layer of Ni is deposited by PVD. These wafers are then heated in an H₂ atmosphere to reduce them, which provokes the formation of nanoparticles through unwetting of the Si surface. These nanoparticles are then used as catalyst for the growth of vertically aligned carbon nanotubes by PECVD. As the diameter and density of the tubes are directly related to the diameter and density of the nanoparticles, we are interested in getting a clear idea of these parameters before performing the nanotube growth.

Some Ni-covered Si wafers were prepared before the first confinement (substrates labeled as old), some were prepared in September (new substrates). Here, we are interested in seeing whether the age and thickness of the Ni layer plays a role on the nanoparticles size and density. Also, the other parameters to study are the temperature at which the unwetting is performed, as well as the duration of this reaction.

To perform this study, we prepared samples from various substrates at various temperatures and during various times. The substrates are then observed with SEM, and several pictures are taken to increase the statistics. These pictures are then analyzed with ImageJ, as shown on Figure 1.1.

Figure 1.1: Typical SEM image of Ni nanoparticles: from the raw image to particle analysis

In this exercise, we are going to treat the tables obtained from ImageJ: these tables contain the (x,y) positions of the particles as well as their area.

2 Data wrangling

• Load the packages tidyverse, readxl, units, broom, and ggforce. We are going to getting used to work with units in this exercise, which is a very good habit to take in order to avoid many unit conversion problems. The package ggforce allows doing ggplot plots with tibbles containing units. Set the global ggplot2 theme to black and white. Also, make it so that the strip.background (background of the facets titles) is blank, and that the strip.text is bold.

library(tidyverse)
library(readxl)
library(units)
library(ggforce)
library(broom)
theme_set(theme_bw()+
          theme(strip.background = element_blank(),
                strip.text=element_text(face = "bold", size=14)
                ))

• Find all MLxx_xx_cc.csv files in the Data folder and store them in flist.

flist <- list.files(path="Data", pattern="m.csv")

• Read the sample.xlsx file that contain all characteristics of the various samples, such as their temperature, substrate type, time of reaction, and store the result in samples.
  • Modify the tibble samples so that its columns are named “sample”, “T”, “time”, and “substrate”
  • Using separate(), separate the “substrate” column into “sub_thick” and “sub_age” containing the thickness and age of the substrate. Use convert = TRUE to convert the characters to integers if applicable.
  • Give the columns their proper unit when applicable (refer to this example)

samples <- read_excel("Data/sample.xlsx") %>% 
    rename(sample    = "name",
           T         = "Temperature (°C)",
           time      = "Time (minutes)",
           substrate = "Substrate (nm)") %>% 
    separate(substrate, c("sub_thick","sub_age"), convert = TRUE) %>% 
    mutate(sub_thick = set_units(sub_thick,"nm"),
           time      = set_units(time,"min"),
           T         = set_units(T,"degC")
           )

• Create the readfile(filename) function that, given one of the csv files names, will:
  • Determine the unit used in this file: from filename, store the unit into the variable UNIT (that should be the string "um" or "nm"). You can use unlist(strsplit(filename,"_")) to get a vector of the elements of filename separated by a _ character.
  • Read the csv file, then successively:
  • Select only the X, Y and Area columns (some files have more columns) and rename them to lowercases names.
  • Add the file column containing the filename, and then separate it into 3 columns sample, number and unit
  • Remove the unit column
  • Apply the corresponding units to x, y and area. You can apply a unit to a vector x based on a string xx using set_units(x, xx, mode = "standard").
  • Create a column diameter containing the diameter of the particles.
• Test this function on 2 files with 2 different units to check that it gives the expected result.

readfile <- function(filename){
    UNIT <- unlist(strsplit(filename,"_"))
    UNIT <- gsub(".csv","",UNIT[length(UNIT)])
    read_csv(filename) %>% 
        rename_all(tolower) %>% 
        select(x,y,area) %>% 
        mutate(file=gsub("Data/","",filename)) %>% 
        separate(file, c("sample","number","unit"), sep="_") %>% 
        select(-unit) %>% 
        mutate(x        = set_units(x, UNIT, mode = "standard"),
               y        = set_units(y, UNIT, mode = "standard"),
               area     = set_units(area, paste0(UNIT,"*",UNIT), mode = "standard"),
               diameter = sqrt(4*area/pi))
}
readfile("Data/ML16_02_um.csv")

## # A tibble: 100 x 6
##        x     y   area sample number   diameter
##     [um]  [um] [um^2] <chr>  <chr>        [um]
##  1 0.221 0.116  0.006 ML16   02     0.08740387
##  2 0.328 0.273  0.009 ML16   02     0.10704745
##  3 0.497 0.423  0.007 ML16   02     0.09440697
##  4 0.189 0.638  0.011 ML16   02     0.11834541
##  5 0.421 1.067  0.010 ML16   02     0.11283792
##  6 0.222 1.410  0.007 ML16   02     0.09440697
##  7 0.634 1.438  0.011 ML16   02     0.11834541
##  8 0.634 1.438  0.011 ML16   02     0.11834541
##  9 0.408 1.576  0.009 ML16   02     0.10704745
## 10 1.677 0.294  0.009 ML16   02     0.10704745
## # … with 90 more rows

readfile("Data/ML16_04_nm.csv")

## # A tibble: 95 x 6
##          x        y      area sample number  diameter
##       [nm]     [nm]    [nm^2] <chr>  <chr>       [nm]
##  1 485.145  176.363  8177.769 ML16   04     102.04048
##  2 552.574  421.864  9433.278 ML16   04     109.59390
##  3 907.114   82.537 14149.916 ML16   04     134.22456
##  4 781.287  193.600  5123.831 ML16   04      80.77044
##  5 762.260  447.889  9433.278 ML16   04     109.59390
##  6 521.046  927.668 16118.010 ML16   04     143.25532
##  7 781.635  653.703  8008.106 ML16   04     100.97642
##  8 906.002 1081.880  9501.143 ML16   04     109.98741
##  9 869.637 1270.324  8381.365 ML16   04     103.30288
## 10 162.661 1767.427 10010.133 ML16   04     112.89507
## # … with 85 more rows

• Using readfile() that you just defined, read all csv files and store them into a tidy tibble called particles. Do not use a for loop to do so. Join this table with the samples one. You will see that, since we attributed units to some columns, all data are automatically converted to a single unit. Filter the data for diameters lower than 40 µm as some very large particles were detected in the image processing that are actually not particles.

particles <- tibble(file=flist) %>% 
    mutate(data = map(file, ~readfile(file.path("Data",.)))) %>% 
    unnest(data) %>% 
    inner_join(samples) %>% 
    filter(diameter<set_units(40,"um"))

In case you didn’t manage to get there, here is the particles tibble (it doesn’t contain the units though as you can’t save it in a text file.)

3 Plotting and analysis

3.1 Size analysis

• Now, plot the histogram of all particle diameters, with a fill color depending on the time (you need to convert time to a factor), and with a grid showing temperature vs. substrate age and thickness. Put the legend on top of the graph, and add some transparency to your colors.
• In fact, I usually prefer to plot it using geom_density() which is basically an histogram convoluted with a Gaussian distribution of bandwidth bw. This allows for smoother graphs. Make this plot and play with the bw parameter.
• Convert – with ggplot – the unit of the particle diameters to nanometers or any other unit you want.

particles %>% 
    ggplot(aes(x=diameter, fill=factor(time)))+
        geom_density(alpha=.5, color=NA, bw=4)+
        labs(x    = "`Particle Diameter`",
             y    = "Density [arb. units]",
             fill = "Time [min]")+
        facet_grid(reorder(paste(T,"K"),T)~
                   reorder(paste(sub_thick,"nm Ni -",sub_age),sub_thick), 
                   scales="free_y")+
        theme(legend.position = "top")+
        scale_x_unit(unit="nm", limits=c(0,300))

• Now, store in particles_ave the average particle diameter and its standard deviation per substrate thickness and age, time and temperature of reaction. You will note that mean() keeps the unit of vectors while sd() loses it. Make sure that the standard deviation column has the proper unit (use units(a) <- units(b)).

particles_ave <- particles %>% 
    group_by(sub_thick, sub_age, time, T) %>% 
    summarise(diam   = mean(diameter),
              sddiam = sd(diameter))
units(particles_ave$sddiam) <- units(particles_ave$diam)

• Plot the average diameter evolution with reaction time, with a color per substrate thickness, and on a grid showing substrate age vs temperature.
  • Don’t forget to add error bars corresponding to the standard error of the diameters distribution.
  • Add a line showing a linear fit for all groups.
  • Make sure both plot axes go to 0.
  • Put the legend on top.
  • Define nice axis labels (with sentences instead of column names). In case a column has a unit that was attributed, it will return an error in case the axis label you want to give contains white spaces. You will thus need to use backticks, like so: "`Two words`"

particles_ave %>% 
    ggplot(aes(x=time, y=diam, color=factor(sub_thick)))+
        geom_point(alpha=.5)+
        expand_limits(x = 0, y=0)+
        geom_errorbar(aes(ymin=diam-sddiam,ymax=diam+sddiam), width=.2)+
        geom_smooth(method="lm", se=FALSE, )+
        facet_grid(T~sub_age)+
        labs(x     = "Time",
             y     = "`Average particle diameter`",
             color = "Substrate thickness [nm]")+
        theme(legend.position = "top")+
        scale_y_unit(unit="nm")

• Using broom display the slopes and intercept of all linear fits, without using a for loop. This doesn’t work well with units, so prior to doing the fit, remove the units of your time and diameter columns using as.vector().

particles_ave %>% 
    mutate(diam   = as.vector(set_units(diam,"nm")),
           sddiam = as.vector(set_units(sddiam,"nm")),
           time   = as.vector(time)
           ) %>% 
    nest(data=-c(sub_thick,sub_age,T)) %>% 
    mutate(fit    = map(data, ~lm(data=., diam~time, weights = 1/sddiam)),
           tidied = map(fit, tidy)) %>% 
    unnest(tidied) %>% 
    select(sub_thick,sub_age,T,term,estimate,std.error) %>% 
    mutate(term=gsub("\\(Intercept\\)","y0",term)) %>% 
    mutate(term=gsub("time","slope",term)) %>% 
    pivot_wider(names_from = term, 
               values_from = c(estimate,std.error)) %>% 
    set_names(~ str_to_lower(.) %>%
                str_replace_all("estimate_","") %>% 
                str_replace_all("std.error_","d")) %>% 
    select(sub_thick,sub_age,t,y0,dy0,slope,dslope)

## # A tibble: 6 x 7
## # Groups:   sub_thick, sub_age [3]
##   sub_thick sub_age     t    y0    dy0  slope  dslope
##        [nm] <chr>    [°C] <dbl>  <dbl>  <dbl>   <dbl>
## 1         5 new       750  80.5   4.36 -1.31    0.270
## 2         5 new       850  96.2 NaN    -1.48  NaN    
## 3         5 old       750 191.   14.1  -5.80    0.878
## 4         5 old       850  92.8 NaN    -0.210 NaN    
## 5        10 old       750 198.   12.4  -3.56    0.775
## 6        10 old       850 159.  NaN    -2.35  NaN

3.2 Density and ordering analysis

• Now we want to see the evolution of the density of particles and their ordering. One way of doing this is to look at the G(r) distribution, i.e. the probability to find a particle in r if one is in 0. I provide here below the function gofr(x,y,dr,Rmax) that computes it between 0 and Rmax with a step dr, provided the (x,y) positions of particles.

gofr <- function(x, y, dr=.2, Rmax=10){
    # Make sure units are uniform before using dist()
    units(y) <- units(x) 
    # dr and Rmax are unitless but should be given in the same units as x
    dr   <- as.vector(dr)
    Rmax <- as.vector(Rmax)
    # Get a vector of all Euclidian distances
    dd   <- as.vector(dist(tibble::tibble(x,y)))
    # Make a histogram out of it
    dd.hist <- hist(dd, 
                    breaks=seq(0, max(dd)+dr, by=dr),
                    plot=FALSE)
    # Get the r values
    r <- dd.hist$mids
    # Compute the normalization by the surface of the 
    # ring of radius r and thickness dr
    rlo  <- r - dr/2
    rup  <- r + dr/2
    ring <- pi*(rup^2 - rlo^2)
    # Return the tibble containing r and G(r) with the same unit as x
    # Only data for r<Rmax is wanted, and we remove the first element too
    d <- tibble::tibble(R    = r[r<Rmax], 
                        GofR = dd.hist$counts[r<Rmax]/ring[r<Rmax]) %>% 
            rename(r="R",gofr="GofR")
    units(d$r) <- units(x)
    d[-1,]
}

• Using pipe operations, compute G(r) for each image and store it into particles_gofr.
  • We want to compute G(r) up to 500 nm by step of 2 nm. Make sure x, y, dr and Rmax are given with the same unit.
  • Then, compute the average G(r) for each sample.
  • Join particles_gofr with samples to retrieve the samples information

particles_gofr <- particles %>% 
    mutate(x = set_units(x,"nm"),
           y = set_units(y,"nm")
           ) %>% 
    nest(data = -c(sample, number)) %>% 
    mutate(GofR = map(data, ~gofr(.$x, .$y, dr=2, Rmax=500))) %>% 
    unnest(GofR) %>% 
    group_by(sample, r) %>% 
    summarize(gofr=mean(gofr)) %>% 
    inner_join(samples)

• Find the best way to represent G(r) for all samples, that allows seeing the evolution with all parameters.

particles_gofr %>% 
    ggplot(aes(x=r, 
               y=gofr, 
               color=factor(time)))+
        geom_smooth(method="loess", alpha=.5, span=.15, se=FALSE)+
        geom_line(alpha=.2)+
        expand_limits(x = 0, y=0)+
        facet_grid(reorder(paste(T,"K"),T)~
                   reorder(paste(sub_thick,"nm Ni -",sub_age),sub_thick), 
                   scales="free_y")+
        theme(legend.position = "top")+
        labs(x     = "r",
             y     = "G(r)",
             color = "Time [min]")